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ABSTRACT 


This  paper  presents  a  new  set  of  differential  methods  for  solving  the  inverse 
scattering  problem  associated  to  the  propagation  of  waves  in  an  inhomogeneous 
medium.  By  writing  the  medium  equations  in  the  form  of  a  two-component  sys¬ 
tem  describing  the  interaction  of  rightward  and  leftward  propagating  waves,  the 
causality  of  the  propagation  phenomena  is  exploited  in  order  to  identify  the 
medium  layer  by  layer.  The  recursive  procedure  that  we  obtain  constitutes  a 
continuous  version  of  an  algorithm  first  derived  by  Schur  in  order  to  test  for  the 
boundedness  of  functions  analjd.ic  inside  the  unit  circle.  It  recovers  the  local 
reflection  coefficient  function  of  the  medium.  Using  similar  ideas,  some  other 
differential  methods  are  also  derived  to  reconstruct  alternative  parametriza- 
tions  of  the  layered  medium  in  terms  of  the  local  impedance  or  of  the  potential 
function.  One  of  these  methods  is  known  in  the  literature  as  the  method  of 
characteristics. 

The  differential  inverse  scattering  methods  turn  out  to  be  very  efficient 
since,  in  some  sense,  they  let  the  medium  perform  the  inversion  by  itself  and 
thus  fully  exploit  its  structure.  They  provide  an  alternative  to  classical  methods 
based  on  integral  equations,  which,  in  order  to  exploit  the  structure  of  the  prob¬ 
lem,  must  ultimately  resort  to  differential  equations  of  the  same  type. 


1.  INTRODUCTION 


The  inverse  problem  for  the  one-dimensional  Schrodinger  equation  and  for 
two-component  scattering  systems  has  received  a  large  amount  of  attention 
over  the  years.  This  interest  is  motivated  by  the  numerous  applications  of  such 
problems  existing  in  fields  as  varied  as  geophysics,  transmission-line  analysis, 
filter  design,  voice  S3mthesis  and  quantum  physics  [l]-[lO]. 

The  first  complete  solution  of  the  inverse  scattering  problem  was  obtained 
by  Gelfand  and  Levitan  [ll],  in  the  context  of  reconstructing  a  second  order 
differential  operator  from  its  spectral  function.  Subsequently,  several  alterna¬ 
tive  solutions  were  proposed  by  Marchenko  [12],  Krein  [13],  Kay  and  Moses  [14] 
and  Faddeev  [8], [15].  Other  inversion  procedures  were  derived  by  Gopinath  and 
Sondhi  [5], [6]  and  by  Zakharov  and  Shabat  [16],[17]  for  systems  described 
respectively  by  transmission-line  type  equations  and  by  two-component  scatter¬ 
ing  models. 

Since  all  the  Inverse  scattering  procedures  mentioned  above  were  formu¬ 
lated  in  terms  of  integral  equations  it  was  widely  accepted  in  the  scientific  com¬ 
munity  that  inverse  problems  require  the  solution  of  such  equations.  However, 
independently  of  the  work  of  mathematicians  and  physicists,  geophysicists  such 
as  Goupillaud,  Claerbout  and  Robinson  developed  approaches  v.'hich  more 
directly  exploit  the  physical  properties  of  layered  media  in  w’hich  waves  pro¬ 
pagate.  Their  solutions,  sometimes  referred  to  as  dynamic  deconvolution 
methods  [2], [4],  reconstruct  the  medium  layer  by  layer,  in  a  recursive  manner. 
However  this  work  was  formulated  in  terms  of  a  discretized  layered  earth  model 
and  was  therefore  not  recognized  as  providing  a  solution  to  the  general  immerse 
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scattering  problem,  In  fact,  when  dealing  with  continuously  varying  media  geo¬ 
physicists  went  back  to  using  integral  equations  based  approaches  [l],[3],[l8]. 

More  recently,  Deift  and  Trubowtz  [19]  proposed  a  potential  reconstruction 
method  based  on  a  trace  formula  which  calls  for  the  propagation  of  an.  ordinary 
difierential  equation  and  which  does  not  fit  the  classical  inverse  scattering 
framework. 

The  objective  of  this  paper  is  to  give  a  comprehensive  account  of 
difierential  inverse  scattering  methods.  This  is  done  by  first  deriving  an 
infinitesimal  layer  peeling  procedure  which  can  be  viewed  as  a  continuous  ver¬ 
sion  of  the  dvmamic  deconvolution  algorithm.  This  algorithm  is  in  fact  a  continu¬ 
ous  form  of  a  method  used  by  Schur  to  test  for  the  boundedness  of  functions 
analytic  inside  the  unit  circle  [20], [2i].  The  identification  of  recursive  layer 
extraction  methods  -with  the  Schur  algorithm  was  first  made  by  Dewilde  and  his 
coworkers  [22], [23].  The  method  of  characteristics  (see  e.g.  Sym.es  [24]),  used 
by  Santosa  and  Schwetlick  [25]  and  by  Sondhi  and  Resnick  [26]  for  sohdng 
acoustical  inverse  problems,  can  also  be  interpreted  from  this  point  of  view.  Tne 
relation  between  the  difierential  inverse  scattering  methods  that  we  propose 
and  the  classical  integral-equations-based  approaches  is  then  discussed.  It  is 
shown  that  by  exploiting  the  structure  of  these  integral  equations,  one  can 
obtain  a  system  of  difierential  equations  which  solves  the  inverse  problem.  The 
difierential  equations  have  the  same  dynamics  as  the  Schur  recursions  but 
reo^uire  certain  boundary  values  that  have  to  be  successively  computed  by  using 
the  integral  equations.  In  fact  these  recursions  are  of  the  same  type  as  the 
Kre  in-Levins  on  equations  for  factoring  the  resolvent  of  a  Toeplitz  operator 
[27],[2B]. 

The  paper  is  organized  as  follows.  Several  physical  models  of  a  scattering 
medium  are  presented  in  Section  1.  These  provide  various. equivalent  parametr- 


izations  of  the  medium  and  give  rise  to  different  formulations  of  inverse  scatter¬ 
ing  problems.  The  continuous  layer-peeling  algorithm  and  the  associated  Schur 
recursions  for  reconstructing  the  reflectivity  function  parametrization  of  the 
medium  are  derived  in  Section  3.  They  are  then  used  in  Section  4  to  obtain 
other  differential  methods  that  reconstruct  either  the  local  impedance  or  the 
Schrodinger  potential  medium  parametrizations.  Section  5  relates  these 
differential  methods  to  the  integral  equations  approaches  and  describes  the 
Krein-Levins on  ty'pe  differential  solution  of  the  inverse  problem.  In  Section  6  the 
results  of  the  earlier  sections  are  extended  to  some  cases  when  the  scattering 
medium  is  not  lossless  and  Section  7  concludes  with  observations  on  possible 
extensions  of  these  results. 
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2.  Physical  Models  of  Scattering  Media 

The  inverse  scattering  methods  that  we  discuss  in  this  paper  concern 
several  classes  of  physical  models  which  correspond  to  equivalent  descriptions 
of  a  lossless  scattering  medium.  They  arise  in  the  study  of  transmission-lines 
and  of  vibrating  strirrgs,  in  the  analysis  of  layered  acoustic  media  and  of  the 
vocal  tract  and  in  the  description  of  particle  scattering  in  quantum  physics  [l]- 
[9],[29],[303. 

The  first  model  that  we  consider  is  described  by  the  symmetrized 
telegrapher's  equations 

u{x.t)  0  —Z{x)-^  v{x,t) 

which  may  be  viewed  as  obtained  from  the  usual  transmission-line  equations  by 
assuming  that  the  inductance  per  unit  length  equals  the  inverse  of  the  capaci¬ 
tance.  Z{x)  in  the  above  equation  corresponds  to  the  local  impedance  for  a 
transmission-line  or  to  the  area  function  of  the  vocal  tract  model  [6], [25], [26]. 
Since  in  equation  (2.1)  the  "voltage"  and  "current"  variables  are  expressed  in 
different  units,  we  also  consider  the  normalized  quantities 

V{x  ,t)  =  v{x  ,t)Z{x)~^''^  and  I{x  ,t)  =  i{x  ,t)Z{x)'^''^  (2.2) 

which  now  have  the  same  dimension.  In  terms  of  these  normalized  variables 
(2.1)  becomes 


6 

V{x,t)' 

dx 

-A:(s) 


0 

dt 


dt 
k  [x) 


V{x,t) 


(2.3) 


where  A:(ar)  is  the  local  reflection  coefficient  (also  called  the  reflectiznly  June- 
tion)  given  by 


k{z)  =  Z{x)  Z(x) 


^  _  .1  ^ 


2  dx 


—In  Z{x) 


(2.4) 


Note  that,  as  a  direct  consequence  of  this  normalization,  we  have 

i{x.t)  I{x.t)  ^  • 


(2.5) 


From  the  system  (2,3)  we  can  obtain  directly  the  second  order  -wave  equations 


(^-^)/(^.0  -  9(:^)  /(^.O  =  0 


(2.6) 


where  the  potentials  are  given  by 

d 


P{x)  =  —■—^{x)+k{x)'^  =  Z(x)i/2 — 1/2 
dx  dx'^ 


Q(z)  =  ^(z)  +  kixr  =  ^2(x)./s 


(2.7) 


In  the  transform  domain  the  equations  (2.6)  take  the  form  of  Schrodinger  equa¬ 
tions.  which  justifies  calling  P{x)  and  Q{x)  potentials.  From  (2.3)  we  can  also 
obtain  a  model  where  the  variables  of  interest  are  right  and  left  propagating 
waves  defined  as 

Wjiix.t)  -  ^  and  Wiix.t)  =  g) 


The  evolution  of  the  wave  variables  is  given  by 


-  B-- 


6 

dx 

6 

dt 

-k  (z) 


-k{x) 

a 

dt 


WpXx,t) 

VfL{x,t) 


(2-S) 


To  iaterpret  this  equatioa,  note  that  when  the  impedance  is  constant  over  a  cer¬ 
tain  section  of  the  medium  we  shall  have  k{x)  =  0  and  therefore 
Wp{x,t)x  Wp{t—x)  and  Wi{x,t)-  Wi{t-^x),  corresponding  to  non-interacting 
right  and  left  propagating  waves.  The  intensity  of  local  interaction  between  the 
waves  propagating  in  opposite  directions  is  quantified  by  A:(z),  which  justifies 
calling  it  the  local  reflection  coefficient.  A  simple  discretization  of  (2.9)  gives  the 
lattice  model  sho'^vn  in  Fig.  1.  Such  discrete  lattice  structures  appear  in  a  large 
number  of  applications,  such  as  the  linear  prediction  algorithms  for  speech  sig¬ 
nals  [31],  the  layered-earth  models  of  Goupillaud  [2], [3],  and  digital  filter  syn¬ 
thesis  [22].  The  model  in  Fig.  1  is  in  fact  crucial  to  the  intuitive  understanding 
of  the  inverse  scattering  techniques  that  we  shall  derive  below. 

By  performing  the  space  transformation 

y{x)  =  fZiOdi  (2.10) 

0 

on  the  telegrapher’s  equation  we  obtain  the  string  eqvxition 

■^v{y,t)  =  t^{y)^v{y  ,t)  (2.11) 

where  juiy)  is  the  mass  density  of  the  string  and  is  given  by 

fi{y)  =  Z-^[x(y)]  (2.12) 

where  x{y)  is  the  inverse  transformation  corresponding  to  (2.10).  This  model 
arises  in  connection  with  the  use  of  inverse  scattering  methods  in  linear  estima¬ 
tion  theory  [29],  By  using  the  alternate  space  transformation 
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y'{z)  =  JZiiVdi  (2.13) 

0 

we  also  obtain  the  conjugate  string  equation 

-^i{y\t)  =  ^'{y')-^iiy\t)  (2.14) 

with 

fx'iy')  =  Z^lz(y')]  ■  (2-15) 

Note  that  /zly(x)]  ju'ly'(x)]  =  1,  which  explains  referring  to  (2.11)  and  (2.14)  as 
conjugate  equations. 

The  four  models  of  a  scattering  medium  that  we  use  in  the  sequel  are  thus 
the  telegrapher's  equations  (2.1)  parametrized  by  Z(x),  the  Schrodinger  equa¬ 
tions  (2.6)  parametrized  by  P(3:)and9(x).  the  two-component  wave  system  (2.9) 
specified  by  A:(z)  and  the  string  equations  (2.11)  and  (2.14)  parametrized  by 
/ii(y)  and  fZ[y')  respectively.  The  objective  of  the  inverse  scattering  problem 
that  we  address  below  is  to  reconstruct  any  of  the  above  parametrizations  from 
some  given  scattering  data.  The  scattering  data  is  obtained  by  probing  the 
medium  in  order  to  determine  its  impulse  or  frequency  response  at  one  of  the 
boundaries.  The  probing  signals  and  the  medium  response  are  assumed  to  be 
measured  perfectly,  i.e.  the  scattering  datauhll  be  considered  noise  free.  Also 
note  that,  since  k{x)  and  P(x)and9('2:)  are  expressed  in  terms  of  the  first  and 
second  derivatives  of  the  impedance  function,  the  different  inversion  methods 
will  require  various  degrees  of  smoothness  for  Z[x).  Yfhen  discussing  the  vari¬ 
ous  cases  we  therefore  assume  that  the  Z{x)  function  is  as  smooth  as  necessary 
for  the  expressions  involved  to  be  well-defined;  fairly  standard  limiting  pro¬ 
cedures  can  often  be  used  to  relax  these  restrictions. 

The  inverse  scattering  problem  associated  with  the  Schrodinger  equation  of 
quantum  physics  is  complicated  by  the  possible  existence  of  bound  states.  A 
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consequen.ce  of  the  assumed  transmission-line  model  (2.1)  is  that  energy  cannot 
be  trapped  in  the  medium  thereby  ruling  out  the  possibility  of  bound  states 


T 


3.  Continuous  Parameter  Schur  Recursions 


The  basic  differential  inverse  scattering  method  that  we  discuss  in  this 
paper  relies  on  the  wave  picture  associated  ■with  equation  (2.9),  a  discrete 
approximation  of  which  is  depicted  in  Fig.  1. 

3. 1.  The  scattering  data 

The  necessary  data  for  the  reconstruction  of  the  scattering  medium  param¬ 
eters  may  be  obtained  in  two  possible  ways. 

In  the  first  case  the  medium  is  assumed  to  be  quiescent  at  t=0  and  it  is 
probed  by  a  known  rightward  propagating  waveform  incident  on  the  medium 
after  t=0.  This  waveform  F/?(0,i)  ■will  in  general  be  an  impulse  followed  (in  time) 
by  a  piecewise  continuous  function,  but  we  also  discuss  the  case  when  no  leading 
impulse  is  present.  The  measured  data  is  the  leftward  propagating  wave,  as  it  is 
recorded  at  a:=0,  Wi{0,t).  It  can  be  viewed  as  obteiined  by  convolving  the 
impulse  response  R{t).  of  the  scattering  medium,  with  the  probing  wave 
Since  the  ultimate  objective  is  to  measure  the  impulse  response  of  the 
medium,  the  nature  of  the  probing  wave  is  not  important  provided  it  contains 
enough  energy  at  all  frequencies.  Rote,  indeed,  that  as  long  as  is  given 

and  f>'i(0,O  is  measured  perfectly,  we  can  always  obtain  the  impulse  response 
by  performing  a  deconvolution. 

Another  way  of  gathering  scattering  data  is  to  perform  a  measurement  of 
its  frequency  response  Bio)  by  sending  into  the  medium  sinusoidal  waveforms 
at  various  frequencies  and  measuring  the  magnitude  and  phase-shift  of  the 
returning  sinusoidal  wave.  This  is  clearly  equivalent  to  the  time-domain  meas- 
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urements  described  above,  since  R{o)  is  the  Fourier  transform  of  E{t). 

From  a  practical  point  of  view  we  cannot  ahrays  directly  generate  the 
waveform  and  measure  However  we  usually  do  have  access  to 

the  physical  variables  v(D,t)  and  i(O.f),  and  by  obtaining  the  medium  response 
in  terms  of  these  variables  we  can  reconstruct  the  corresponding  and 

fVl(0,t)  by  using  (2.8).  (In  the  sequel  we  assume  that  Z(0)  =  1.)  The  nature  of 
the  measurements  (impulse  or  frequency  response)  clearly  depends  on  the  phy¬ 
sical  apparatus  that  is  available.  In  the  geophysical  context,  approximate 
impulse  responses  are  obtained  by  using  explosive  sources  (dynamite,  air-guns) 
and  frequency  response  data  can  be  generated  by  using  wide-band  acoustic 
sources  [2], [32] 


3.2.  The  layer  peeling  procedure 

Suppose  that  the  incoming  wave  JVj}(0,t)  contains  a  leading  impulse.  This 
impulse  will  propagate  through  the  medium  and,  since  the  medium  is  causal,  it 
is  not  hard  to  recognize  by  examining  Fig.  1.  that  the  waves  }Vz^(x,t)  and  Wi{x,t) 
must  have  the  form 


W}{{x,t)  =  6{t —x)  +  wjf>{x ,t)  u{t —x) 
=  wi{x  ,t)u{t-x) 


(3.1) 


where  \uj^{x  ,t)  and  -wiix  .t)  are  some  piecewise  continuous  functions,  d(-) 
denotes  the  Dirac  distribution,  andu(-)  is  the  unit  step  function,  i.e., 

^  ^  f  1  for  t>0 

0  for  t<0  .  (3.2) 


The  causal  nature  of  Wp/yZ,t)  and  Wi{z,t),  i.e.  that  they  are  zero  for  i  "<  z.  is  a 
direct  consequence  of  the  fact  that  the  medium  was  at  rest  at  f  =0  since  the 
impulse  requires  an  amount  of  time  equal  to  z  to  reach  the  depth  z  in  the 
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mediuin.  Note  that  we  ass'omed  that  the  perturbation  in  the  medium  originated 
from  its  left  end  alone.  By  substituting  (3.1)  into  the  propagating  equations 
(2.9)  and  equating  the  coefficients  of  ,6(t  -x)  on  both  sides,  we  find  that 

vji(x,x+)  -  |*(z)  and  ~^j,(x,z+)  =  -  ^(x)^  .  (3.3) 


This  argument  is  an  application  of  the  classical  method  of  propagation  of  singu¬ 
larities  (see.  e.g.  [33]).  Now,  noting  that  the  expression  (3.1)  for  Wf>{x .t)  and 
Wi{x,t)  implie s  that 


V{x,t)  =  5{t —x)  + '^{x  ,t)  u{t —x) 
I(x,t)  =  5(f -z)  4- 'l(z ,£)  u(f -z) 


(3.4) 


where 

'ir{x  ,t)  =  wji{x  ,t)  +  u>i{x  ,t)  and  ^{x  ,t)  =  wj>{x  ,t)  -  u>i{x  ,t)  ,  (3.5) 


we  conclude  from  (3.3)  and  (2.7)  that  the  potentials  are  given  by 

P(X)  = 

Q{x)  =  —2~^{x,x+) 


(3.6) 


The  above  results  show  that  the  local  reflection  coefficient  sequence  fc(z) 
can  be  reconstructed  directly  from  the  reflected  waves  at  depth  x  in  the 
scattering  medium.  However  we  have  assumed  that  only  the  reflected  wave  at 
x  =  0  is  measured;  the  waves  at  depth  z  >  0  will  be  constructed  by  a  recursive 
procedvire.  Thus  let  us  assume  that  the  waves  at  point  z  have  already  been 
computed;  then  k{x)  can  be  readily  identified  as  Ti;jr,(z,z+)  using  (3.3)  and  sub¬ 
stituting  this  value  into  the  propagation  equations  (2.9),  we  can  compute  the 
wmves  at  depth  z+A.  Therefore,  starting  at  z=0,  the  resulting  recursive  algo¬ 
rithm  can  successively  identify  the  local  reflection  coefficient  for  increasing 
values  of  z.  This  recursive  inverse  scattering  process  may  also  be  viewed  as  a 
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layer-reeling  algorithm,  where  at  every  step  one  infinitesimal  layer  of  the 
scattering  medium  is  identified  and  effectively  removed.  The  right  and  left  pro¬ 
pagating  waves  inside  the  medium  are  recursively  generated  and  can  be 
regarded  at  each  step  as  a  new  set  of  scattering  data  for  the  remaining  extent  of 
the  medium.  For  a  lossless  and  discrete  layered  medium  this  algorithm  is  known 
in  geophysics  as  a  dynamic  deconvolution  process  [4]  and  it  is  called  the  do^Lfn.- 
xvard  continuation  method  by  Bube  and  Burridge  [34].  Dewilde  et  al.  [22], [23], 
noted  that  this  algorithm  is  equivalent  to  the  Darlington  synthesis  procedure  for 
scattering  functions  and  pointed  out  its  similarity  to  a  result  of  Schur  (1917) 
that  will  be  discussed  in  the  next  section.  In  the  context  of  fast  algorithms  for 
linear  estimation  and  operator  factorization  theory  these  recursions  are  some¬ 
times  referred  to  as  the  fast  Cholesky  recursions.  The  operator  factorization 
identity  associated  to  these  recursions  is  discussed  in  [28]. 

To  implement  the  layer  peeling  method  we  can  use  the  follo^^ing  numerical 
scheme  which  was  also  derived  in  a  slightly  different  form  by  Dewdlde,  Fokkema 
andWidya  [23],  Denote 

=  wy{x  ,t+x)  and  ai[x  ,t)  =  vui{x  ,t  +x)  .  (3.7) 

Then,  integrating  the  evolution  equations  (2.9),  w'e  obtain,  after  some  calcula¬ 
tion,  the  following  system  of  equations 

Z 

aj>{x.t)  =  vJj,{0.t)  -  fk{Oai^{^,t)d^ 

0 

(3.6) 

ai{x,t)  =  wi{0.t -i-Zx)  -  J'ktf)aj^{^.t +2x —Z^)d^ 

0 


together  with  the  formula  giving  the  reflection  coefficient 
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k{z)  =  2aiiz,0)  =  2[vji{0,2z)  -  flcU)aE{i.Zx-Z^)di]  .  (3.9) 

0 

By  recursively  integrating  (3.6)*(3.9)  along  successive  antidiagonals  in  the  {z,t) 
plane,  as  depicted  in  Fig.  2.,  we  can  obtain  the  local  reflection  coefficients  A:(x), 
for  increasing  values  of  z.  Note  also  that  we  only  need  to  know  the  probing  and 
reflected  waves  Wj>{Q,t)  and  over  the  time  span  [0,2x]  in  order  to 

recover  the  transmission  line  parametrization  up  to  depth  z. 

From  a  computational  point  of  view,  if  we  assume  that  the  part  of  the 
medium  of  interest  has  total  length  L.  and  if  we  use  a  difference  scheme  with 
step-size  L/ N  in  the  propagation  of  the  layer  peeling  algorithm,  (3.8)-(3.9),  the 
total  number  of  operations  required  to  reconstruct  the  local  reflection 
coefficient  parametrization  is  0[N^).  These  algorithms  are  therefore  very 
efficient,  when  compared  to  the  direct  deconvolution  methods  w’^hich  do  not 
exploit  the  physical  structure  of  the  medium. 


3.3.  The  Schur  recursions 

By  taking  the  Fourier  transform  of  the  waves  Wj>{z,t)  and  Wi{z,t)  the  pro¬ 
pagation  equations  (2.9)  become 


d 

Wjiix.u) 

-jz 

-k{x) 

dz 

Wi(x.o) 

.-k{x) 

0^ 

(3.10) 


and  the  frequency  response,  or  r2flectiQn  coefficient  function  of  the  section  of 
scattering  medium  over  [a:,'®)  is  given  by  the  ratio 


R{x  ,o) 


(3.11) 


Clearly 
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H{0.u)  =  E{u)  (3,12) 

is  provided  by  the  given  scattering  data.  Using  these  definitions  the  layer¬ 
peeling  algorithm  described  above  can  be  recast  as  a  recursive  procedure  for 
computing  the  sequence  of  reflection  coefficient  functions  R{x,u)  for  increasing 
values  of  x.  Since  R{z,cS)  is  the  ratio  of  variables  with  a  linear  evolution  given  by 
(3.10),  it  will  not  be  surprising  to  find,  after  some  algebra,  that  it  satisfies  the 
Riccati  equation 

-^R{x,cS)  ~  Zj  uR{z  ,cS)  i^{^)[Rix  ,cS)^  -  l]  .  (3.13) 

It  is  not  clear  how  this  can  help,  since  fc(-)  is  unknown,  but  recalling  the  identity 
(3.3)  for  h{x)  and  the  form  (3.1)  for  the  waves  at  x.  we  find  by  using  the  initial 
value  theorem  for  unilateral  transforms  that 

k{x)  =  2u>i{x,z-¥)  =  lim8;'c,>.ff  (x,w)  .  (3.14) 

(j-tm 

By  using  (3.14)  the  equation  (3.13)  can  now  be  propagated  autonomously.  In 
terms  of  the  causal  impulse  response  R{x,t)  corresponding  to  the  reflection 
function  R{x,d),  the  equation  (3.14)  simply  states  that 

kix)  =  ZR{x,Q+)  .  (3.15) 

The  Riccati  equation  (3.13)  for  the  reflection  function  is  fairly  well-known  in  radi¬ 
ative  transfer  and  transmission-line  theory,  and  is  a  direct  consequence  of  the 
rules  of  cascading  infinitesimal  scattering  layers  [35].  More  details  about  the 
evolution  of  the  medium  representation  under  successive  compositions  of 
infinitesimal  scattering  layers  will  be  given  in  Section  5, 

In  the  context  of  the  inverse  problem  of  geophysics  the  Riccati  equation 
(3.13)  "was  also  obtained  by  Gjevik  et  al.  [36].  However,  they  did  not  notice  the 
relation  (3.14)  which  can  be  used  to  propagate  the  Riccati  equation  recursively, 
starting  from  the  scattering  data  (0,£o)  =  R  {u).  They  proposed  an  iterative  , 


rather  than  recursive  procedure  to  compute  the  k{z)  function.  Vfe  should  note 
at  this  point  that  the  computational  issues  associated  with  solving  (3, 13), (3. 14) 
have  not  been'studied  and  deserve  further  investigation. 

The  equations  (3.13)  and  (3.14)  constitute  the  continuous  version  of  a  pro¬ 
cedure  derived  by  Schur,  in  1917  [20], [El],  for  testing  boundedness  of  an  ana¬ 
lytic  function  outside  the  unit  circle  of  the  complex  plane.  Given  a  power  series 
in  in  ^(z),  Schur  proved  that  |S'(z)|  <1  on  the  unit  circle  if  and  only  if  the 
sequence  of  coefficients  generated  by  the  recursion 


^n+I 


1  Sji  ( Z  )  fcyt 

2  i  -_A:„5„(2) 


vAth 


Urn  5„(2) 


(3-16) 


are  such  that  \kn  |  ^  1.  The  discrete  parameter  recursion  (3.16)  is  in  fact  a 
discretized  form  of  the  Riccati  recursion  (3.13)  and  can  be  obtained  from  it  by 
using  a  backwards  difference  scheme. 

The  Schur  algorithm  (3.16)  may  be  interpreted  as  testing  for  the  existence 
of  a  discrete  (i.e.  with  piecewise  constant  impedance)  transmission-line  having 
S'Cz)  for  the  left  reflection  coefficient  function.  Similarly,  the  continuous  version 
of  this  algorithm  may  be  considered  as  testing  for  the  existence  of  a  lossless 
transmission-line  which  synthesizes  the  given  scattering  function  A  con¬ 

dition  for  the  existence  of  such  a  transmission-line  is  that  the  reconstructed 
local  impedance  function  Z{x).  appearing  in  the  model  (2.1),  should  be  strictly 
positive  and  bounded.  Since,  from  (2.4) 


Z{x)  =  .ZiO)  exp\fki^)d^l 
0 


(3.17) 


Z 

this  implies  that  we  need  to  have  \f  A:(|)ci||  <  «=  for  all  x.  In  this  case  E['J)  is 

c 

bounded  by  1  on  the  real  axis.  Vfe  note  that  if  a  transmission-line  is  lossless,  its 
left  reflection  function  B{cj)  must  be  bounded  by  one  on  the  real  axis  as  a  result 
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of  energy  conservation  [5],['7],[9],[l5].  The  above  result  is  therefore  the  con¬ 
tinuous  equivalent  of  the  boundedness  test  devised  by  Schur. 


i  o  - 


4.  Otter  Differential  Inversion  Methods 

In.  the  previous  section  our  analysis  concentrated  on  the  two-component 
system  of  wave  equations  (2.9),  and  in  this  framework  we  have  shown  how^  to 
reconstruct  the  local  reflection  coefficient  function  k  (•).  By  recalling  the  identi¬ 
ties  (2.7)  and  (3.17),  the  potentials  P{  )  and  Q{-),  and  the  impedance  Z{  )  may 
also  be  obtained.  However,  since  fc(-)  is  expressed  as  a  function  of  the  first 
derivative  of  the  local  impedance  function,  the  reconstruction  method  that  we 
have  described  above  requires  the  differentiability  of  Z{  ).  Yfhen  the  local 
impedance  function  is  only  piecewise  differentiable  the  Schur  algorithm  can  be 
modified  to  take  the  discontinuities  into  account.  However  a  more  direct 
method  is  to  use  the  method  of  characteristics ,  which  can  be  described  as  fol¬ 
lows. 


4. 1.  The  method  of  characteristics 

Assume  that  the  probing  wave  Wj>{x,t)  does  not  contain  a  leading  impulse 
and  is  a  piecewise  continuous  function  starting  at  f=0.  Then,  by  causality. 
W^{x,t)  and  Wi{x,t)  must  have  the  form 


Wfi{x,t)  =  wj>{x  ,t)  n{t -x) 
]Vi{x,t)  =  wiix.t)  uit-x) 


(4.1) 


Substituting  these  expressions  into  (2.9)  we  find  that 

wi{x,x+)  =  0  (4.2) 


which  implies  that 
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V{z  ,2+) 


i  /  ON 


Recalling  the  identity  (2.5)  this  sho^vs  that  v;e  have 


V  jx  ,2  +) 

i{x  ,2  +) 


Z{x) 


(4.4) 


Therefore,  to  reconstruct  the  impedance  function,  Z{  ),  we  only  need  to  meas¬ 
ure  the  voltage  and  current  variables  v[0,t)  and  1(0.0  the  left  boundary  of 
the  scattering  medium  and  to  propagate  v(2,0  and  i{x  ,t)  by  using  (4.4)  and 
(2.1)  Note  that  the  knowledge  of  the  voltage  and  current  variables  at  depth  x 
enables  us  to  compute  the  impedance  Z{x),  which  in  turn  can  be  used  to  obtain 
the  functions  v[x-rA,t)  and  i(2+A.O.  In  this  manner  the  impedance  Z[x)  is 
computed  recursively,  starting  from  2  =0  ,  and  this  procedure  is  known  in  the 
literature  as  the  method  of  characteristics  [24]-[26]. 


The  above  inverse  scattering  procedure  can  be  interpreted  in  terms  of  the 
layer-peelirg  technique  of  section  3.2  by  considering  the  discretized  version  of 
(2.1)  shown  in  Fig.  3.  This  figure  indicates  that  the  current  and  voltage  variables 
at  point  (n  +  l)h,  where  A  is  the  discretization  step-size,  are  obtained  from  the 
corresponding  variables  at  depth  nA  by  cascading  a  scattering  layer  described 
by  the  matrix 


S 


n 


Zi■n^y^^  1 

Z(nA) 


(4.5) 


with  tim*e  delays  and  the  inverse  of  the  first  scattering  layer.  This  result  can  be 
obtained  by  noting  that 


lf^(nA,0 

i  (n  A,  f  ) 

=  ©n 

p  (nA.f ) 

where 
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e 


n 


1 

2 


Z{nb)-^^^ 


+  P-][  P-In  +  P:\-^  (4.7) 


is  the  chain  scattering  or  transmission  matri>:-corresponding  to  the  scattering 
representation  [37].  The  projection  matrices  P+  and  P-  appearing  in  the 
above  formula  are  defined  as  follows 


= 


1  0 
P  0, 


P- 


_  fo  o' 
"  lo  i, 


(4.B) 


The  method  of  characteristics  exploits  the  delay  structure  as  displayed  in 
Fig.  3  cind  the  fact  that  the  left  reflection  coefficient  (i.e.  the  2i  entry)  of  the 
matrix  is  the  local  impedance  Z{nA).  Since  both  and  its  inverse  are 
entirely  parametrized  by  the  local  impedance,  the  scattering  layers  associated 
to  these  matrices  can  be  easily  "peeled  off"  (i.e.  their  effect  may  be  accounted 
for)  as  soon  as  Z{nA)  has  been  com.puted. 


The  identity  (4.4)  shows  that  the  reconstruction  procedure  described  above 
can  also  be  used  to  obtain  the  mass  densities  //(•)  and  p'(').  appearing  in  the 
string  equations  (2.11)  and  (2.14).  This  is  done  by  substituting 

12 


piy)  = 


v{y,y+) 


(4.9) 


and 


p'{y')  = 


v(y'  ,y'+) 


i{y'  .1/'+) 


(4.10) 


into  the  equations 


d 

^iy.t) 

o 

1 

v{y,t) 

dy 

/iiy.t) 

^(y.O 

and 


(4.11) 
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a 

V  {x  ,y' ) 

0 

dy' 

>(v'.0 

3 

dt 

0 

^(y',0 

^(y’.O 


that  describe  the  strings  (2.11)  and  (2.14). 


(4.12) 


4.2.  Direct  recovery  of  the  potential 


Similarly,  there  also  exists  a  procedure  for  computing  the  potentials  P{-) 
and  Q{-)  directly,  without  first  reconstructing  the  reflection  coefficient  function 
A:(').  To  do  so  let 


(4.13) 


Then  the  Schrodinger  equations  (8.6)  can  be  rewritten  in  the  form  of  asym¬ 
metric  two-component  differential  systems,  given  by 


_6_|  nx.O]  _ 

0 

dt 

1 

V{x,t) 

P{x) 

a 

dt 

lp(a:,0 

a 

a 

dt 

1 

■  I{x,t)  ' 

Sx  G[x,t) 

Q{^) 

a 

dt 

.  G{x.t)  , 

(4.14) 


(4.15) 


The  layer-peeling  technique  introduced  in  Section  3  can  again  be  used  to 
recover  the  potentials  P(-)  and  Q{-)  directly,  by  noting  that 


Pix)  =  -2F{x,x+) 


(4.16) 


(?(2)  =  —ZG{x,x+) 

Consequently,  if  we  propagate  the  variables  \  V{x  .t),F{x  ,t)l  or  [I{x  ,t),G{x  ,t)l  by 
using  (4.16)  and  the  propagation  equations  (4. 14)-(4.15),.  the  potentials  P(-)  and 


Q(-)  can  be  recovered  directly  from  the  scattering  data. 


To  obtain  the  initial  conditions  for  the  systems  (4.14)  and  (4.15)  we  assume 
that  the  scattering  data  is  W^{0.t)  -  6{t)  and  Wj>{0.t)  =  J?{t)u{t).  Then,  by 
using  equation  (2.3)  and  the  fact  that  k  (0)  =  2I?{0+).  we  find  that 


7(0,0  =  <5(0  + 

F{Q.t)  =  -2[-^(f)  +  ;?(0+)/?(O]  u{t) 


(4.17) 


and 


1(0, t)  =  (5(0  -  i?(0^(0 

G(0,0  =  -2[^(0  +  i?(0+)^(O]^(O 


(4.18) 


"Whereas  in  Section  3  the  potential  was  reconstructed  by  using  the  original 
scattering  data  and  then  differentiating  the  reflectmty  function,  the  method 
that  we  propose  here  first  differentiates  the  scattering  data  and  then  recon¬ 
structs  the  potential  directly. 


The  layer-peeling  algorithm  for  the  systems  (4. 14), (4. 15)  can  be  interpreted 
as  successively  truncating  the  potentials  P{  )  and  Q{  )  in  such  a  ivay  that  the 
new  potentials 


P{z,x)  =  P{z)u{z~x) 
Q{z,x)  =  Q(z)u(z-x) 


(4.19) 


correspond  to  the  part  of  the  original  scattering  medium  located  to  the  right  of 
X.  In  this  interpretation  it  is  assumed  that  the  part  of  the  scattering  medium  on 
the  left  of  X  that  was  removed  by  the  layer-peeling  algorithm  has  been  replaced 
by  free-space  (i.e.  fc(z)  =  0  for  z  <  x).  The  idea  of  using  truncated  potentials  for 
the  analysis  of  direct  scattering  phenomena  was  exploited  earlier  by  Bellman 
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and  Vi'ing  [38]  aaid  is  discussed  in  Lamb  [39].  This  approach  may  also  be 
regarded  as  an  invariant  imbedding  method. 

The  differential  method  presented  above  for  the  reconstruction  of  the 
potentials  P(-)  and  Q{-)  seems  to  be  related  to  the  trace  -method  of  Deift  and 
Trubowitz.  Their  method  is  based  on  the  recursive  computation  of  the  Jost  solu¬ 
tion  of  the  Schrodinger  equation  given  by 

{x,u)  ■¥  [d^  -  P{x)]f  {x,cS)  =  0  (4.20) 

with  boundary  condition 

lim/.(a:,£.>)exp[-y  03:  ]  =  1  .  (4.21) 

Then,  by'  substituting  the  trace  formula 

P{x)  =  —fjuB{u)f^{x.cS)do  (4.22) 

into  (4.20),  f  {x  ,cS)  and  P{x)  can  be  computed  recursively  for  decreasing  values 
of  X .  The  connection  between  the  approach  of  Deift  and  Trubowitz,  and  the  algo¬ 
rithm  that  -we  have  discussed  above  is  not  yet  completely  understood. 
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5.  Integral  Equations  Formulation 

In  Section  3  the  Schur  recursions  were  derived  by  only  using  causality  and 
the  differential  description  of  the  medium.  However  most  classical  inversion 
methods  are  formulated  in  terms  of  integral  equations.  The  objective  of  this  sec¬ 
tion  is  to  derive  a  set  of  Marchenko  integral  equations  for  the  two-component 
system  (2.9)  and  to  show  that  these  equations  can  be  solved  efficiently  by  a  set 
of  differential  equations  similar  in  form  to  the  Schur  recursions. 


5.1.  Transmission  emd  scattering  descriptions  of  the  medium 

The  system  (2.9)  describes  the  transmission  of  waves  through  an 
infinitesimal  section  of  the  medium.  These  infinitesimal  layers  may  be  aggre¬ 
gated  over  the  interval  [0,z]  and  by  using  the  linearity  of  the  medium  we  find 
that  the  waves  Wji[x,t)  and  Wi(z,t)  at  depth  z  are  related  to  the  waves  at  the 
right  boundary  by 


WR{x,t)' 

/i/is(2,f) 

Wdx.t) 

J 

Mziix.t) 

Mzzix.i') 

J 

^i(o.O 

where  *  denotes  the  convolution  operator.  The  matrix 


(5.2) 


is  the  transition  matrix  of  the  medium  over  [O.z]  and  it  satisfies  the  differential 
equation 


M{x.t) 


(5.3) 


with  initial  condition 
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5{t) 

0 


0  ■ 

Sit) 


(5.4) 


The  aggregated  mediiam  corresponding  to  can  be  viewed  as  obtained  by 

composing  the  infinitesimal  layers  that  were  peeled  of!  from  the  medium  by  the 
Schur  algorithm  over  the  interval  [O.z].  Let  M{x,o)  be  the  Fourier  transform  of 
M{x,t).  Then,  the  composition  procedure  for  generating  M  [x /J)  and  the  layer 
peeling  method  are  compared  in  Fig.  4. 


Instead  of  using  the  transmission  description  of  the  medium  given  by  (5.1) 
we  can  use  an  equivalent  scattering  description  which  relates  the  outgoing  waves 
to  the  incoming  waves,  as  follows 


Wp{x,t) 

Tiix.t) 

Pp{x,t) 

WpiO.t)  ' 

WpiO.t) 

Tuix.t) 

The  Fourier  transform  5(z,o)  of  the  matrix 


S{x  ,t) 


Tiiz.t)  Rfiix.t) 
Frix.t)  Tj^ix.t) 


(5.6) 


is  the  scattering  mairix  associated  to  the  medium  over  [0,x]  and  it  can  be 
obtained  from  M[x  ,o)  by  the  relation 

S{x/S)  =  [PJJ{x.o)+P:\[PJ1{x,u)^p;\-\  .  (5.7) 


The  general  rules  of  composition  of  scattering  layers  are  described  in  Redhefler 
[35]. 

As  a  consequence  of  the  delay  structure  and  losslessness  of- the  elementary 
(infinitesimal)  scattering  layers  described  in  Fig.  1,  the  scattering  matrix  S{x,t) 
is  such  that 


and  it  is  lossless,  i.e. 


RpXx.t)  -  Riix.t)  =  0  for  f  <  0 
Tpix.t)  =  Ti{x,t)  =  0  for  t  <x 


(5.6) 

(5.9) 
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S^{z/^)  Six,'^)  =  I 


(5.10) 


^^"here  the  superscript  H  denotes  Hermitian  transpose.  In  the  transmission 
representation  domain  the  relations  (5.6)  and. (5. 9)  imply  that  the  entries  of 
Mix,-)  have  all  support  over  [—x.x].  Finally,  by  noting  that  the  transmission 
medium  is  invariant  when  the  right  and  left  propagating  waves  are  interchanged 
and  time  is  reversed,  we  get  the  following  useful  identities 


Mn{x,t)  =  Mzz{x-t) 
Mzxix.t)  =  M,z{x-t) 


(5.11) 


5.2.  The  MarchenJko  integred  equations 

Vi'hen  the  medium  is  probed  from  the  left,  a  consequence  of  its  delay  struc¬ 
ture  is  that 


WpXz,t)  =  Wi{x,t)  =  0  for  t  <x  .  (5.12) 

By  substituting  (5.12)  into  (5.1)  and  recalling  that  M{x.  )  has  support  on  [-x,x], 
we  obtain  the  system  of  integral  equations 

t  t 

J~  Wp(0,t  -r)Mii{x ,r)dT  +  J" Wi{Q.t  —r)Mizix ,T)dT  =  0 

-s  -r 

(5.13) 

t  t 

f  Wp{Q.t-r)Mzi{x,r)dT  ^  f  -r)Mzz{z .r)dT  =  0 


which  relates  the  entries  of  M{x,t)  to  the  measured  waves  and 

i.e.  the  scattering  data.  From  the  differential  equation  (5.3)-(5.4)  satisfied  by 

M{x,t),  it  can  be  shown  that  Mnix.t)  and  can  be  expressed  as 


Mii{x,t)  =  (5(2-0  ~  ■^(^+0] 

Mzi{x,t)  -  ■mziix,t)[u{x~t)  ~u{x+t)]  . 


Then,  by  using  (5.14)  and  the  symmetry  relations  (o.ll),  we  get  the  Marchenko 
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integral  equations  (for  —x-<t  <r) 

t  r 

J 'Wji{B,t-~T)tnii{x  ,r)d'T  +  f  Wi{0,t +r)Tn2i{x  ,r)dr  =  0 

-X  ~i 

(5.15) 

X  t 

Wi{0.t+s:)  +  f  Wi{0.t+r)mn{x.r)dT  +  f  Wj;{0.t -r)m.zi{x  .r)dr  =  0 

-X 

which  need  to  be  solved  for  the  functions  7nii(x,-)  and  m,gj(2,').  To  guarantee 
the  existence  of  solutions  to  these  equations,  it  is  as  usual  assumed  that  the 
probing  wave  1/^(0. i)  contains  a  leading  impulse,  see  e.g.  (3.1).  In  this  case  the 
integral  equations  (5.15)  take  the  form  of  a  system  of  coupled  Fredholm  equa¬ 
tions  of  the  second  kind 

(  X 

mii{x  ,t)  +  J'wj}{Q,t —r)mii{x  ,r)dT  +  J‘uJi{Q,t+T)mzi{x ,r)dT  =  0 

-z  -i 

(5.16) 

Z  t 

wi(0,t +x)  +  mziix  ,t)  +  J''iJ^i{0,t+r)m.ii{x,T)dr  +  J'w/f{0,t—~)m2i{x,T)dT  =  0 

-t  -X 

The  solution  of  these  equations  can  be  used  to  reconstruct  the  medium  since 
from  (5.3)  and  exploiting  the  form  (5.14)  of  Mii{x,t)  and  Mziix.t)  we  find  that 

*(z)  =  -2m2i(z.x-)  (5.17) 

and 

*®(2)  =  Z-^vn.^^{x,x-)  .  (5.16) 

The  equations  (5.16)  can  be  solved  directly  by  using  a  simple  discretization 
scheme.  If  the  interval  [— 2,z]  is  divided  into  N  equal  subintervals,  this  scheme 
would  require  0{N^)  operations  in  order  to  reconstruct  the  reflection 
coefficient  function  over  [0,z).  •  ■ 

However  the  kernels  iij^(0,f+T)  and  'Wi{B,t-r)  which  appear  in  (5.16)  have 
respectively  a  Toeplitz  and  Hankel  structure  which  can  be  exploited  to  reduce 
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the  number  of  computations.  Thus,  note  that  muiz.t)  and  satisfy  the 

differential  system 


d 

mii(x.t) 

6z 

mzx(z.t) 

_±_ 

dt 

-k{x) 


■m-xiiz.t) 

mzi{x,t) 


(5.19) 


for  —z<t-<x,  -with  initial  conditions 


mjj(0,0)  =  m2i(0,0)  =  0  (5.20) 

IVhen  propagating  (5.19),  it  turns  out  that  we  need  to  supply  the  values  of  the 
kernel  77135(3;, •)  at  t=x—  (providing  k  {x))  and  also  the  value  of  77155(3;, •)  at  t=~x. 
By  using  (5.17)  and  (5.18),  A;  (3;)  can  be  expressed  as 

k{x)  =  77i3i(3;,3;-)  = 

S  Z 

—  j'wi{0.x-¥r)m'^^-^{x  ,r)dr  -  J'WpX0,x-r)m.z\{z ,r)dr  ,  (5.21) 


Fhrthermore  setting  f  =  — 3;  in  (5.16)  we  find  that 

77155(3:, -3;)  =  0  .  (5.22) 

The  differential  system  (5.19),  with  the  boundary  conditions  (5.21)  and  (5.22), 
can  now  be  used  to  compute  77iii(x,f)  and  77135 (3;, f)  recursively.  The  equations 
(5.19)  have  the  same  form  as  the  Schur  recursions  that  were  derived  in  Section 
3.  However  the  Schur  algorithm  is  formulated  as  an  initial  value  problem 
whereas  the  recursions  derived  above  constitute  a  boundary  value  problem. 
These  recursions  are  similar  to  the  Krein-LevinsoTL  equations  for  factoring  the 
resolvent  of  a  Toeplitz  kernel  [27], [2-8].  They  require  the  same  order  of  computa¬ 
tions  as  the  Schur  recursions.  The  stability  of  numerical  schemes  for  propagat¬ 
ing  these  two  types  of  algorithms  is  discussed  in  Gohberg  and  Koltracht  [40]. 

The  differential  equations  (5.19)  could  have  been  derived  also  by  applying 
the  displacement  operators 
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_d_ 

dt 


and 


(5. 23) 


to  the  integral  equations  (5.16),  and  by  using  the  displacement  property 

rf{x-t)  =  0  and  '6-f{x+t)  =  0  (5.24) 

of  Toeplitz  and  Hankel  kernels.  These  displacement  properties  have  been 
exploited  in  [27], [28]  and  [40]  to  derive  some  fast  algorithms  for  the  computa¬ 
tion  of  resolvents  of  Toeplitz  and  Hankel  operators  and  to  obtain  triangular  fac¬ 
torizations  of  these  operators.  Anderson  and  Kailath,  [41],  also  pointed  out  that 
the  resulting  fast  algorithms  can  be  used  to  solve  the  integral  equations  of  clas¬ 
sical  inverse  scattering  theory. 

5.3.  Relation  to  classical  inverse  scattering 

The  integral  equations  (5.16)  are  expressed  in  terms  of  the  general  scatter¬ 
ing  data  Wff{0,t)  =  d{t)  +  and  In  the  litera¬ 

ture,  two  choices  for  the  probing  waves  have  been  made  either  directly  or  impli¬ 
citly 

(i)  wpX^.t)  =  0  and  -  R(t)u{t)  (5.25) 

(ii)  ivpiO.t)  =  wi{0,t)  =  h{t)  (5.26) 

The  second  choice  above  arises  naturally  when  the  scattering  medium  is  ter¬ 
minated  at  its  left  boundary  wth  a  perfect  reflector.  The  Marchenko  integral 
equations  that  we  have  obtained  above  are  therefore  slightly  more  general  than 
those  presented  in  the  literature  of  two-component  inverse  scattering  problems 
[17], [39].  Furthermore  they  can  be  used  to  obtain  the  classical  integral  equa¬ 
tions  that  solve  the  inverse  scattering  problem  for  the  one-dimensional  Schrod- 
inger  equation.  To  do  so,  denote  by 


(5.27) 


w  ^ 

K(x.t)  =  Triiiiz  .t)  -r  m2iix  ,t) 

Thea,  by  adding  the  two  integral  equations  (5.16),  we  obtain 

Z 

VJliO.t +x)  +  K{x  .t)  +  -r)  +  xoi{t —r)]  K{x  ,T)dT  =  0  (5.28) 

-Z 

where  the  potential  is  given  by 

P(x)  =  2  -^(x.x-)  .  (5.29) 

In  the  special  case  when  the  scattering  data  is  given  by  (5.25),  the  above  equa¬ 
tions  correspond  to  the  "classical"  Marchenko  solution  of  the  inverse  scattering 
problem  [9], [15], [33].  When-  (5.26)  is  given  as  scattering  data,  it  can  be  sho'wn 
from  (5.28)  that  the  symmetrized  kernel 

+  Kix.-t)  ]  (5.30) 

satisfies  the  G  elf  and -Levitan  equation  [lO], 

Ks{x,t)  +  |-{/l(2■^-^)  +hix-t)]  + 

*  .. 

f  5-{^(  1 1  “T I )  +  ^(  U  +' I )]  Ks{x ,T)dT  =  0  0  <  f  <  z  (5.31) 

0  * 

and  again 

P{x)  =  2-^(z,z+)  .  (5.32) 

Note  that  the  symmetric  kernel  is  half  the  sum  of  all  the  entries  in  the 

transmission  matrix /^(z.f). 

Finally,  if  we  define 

L{x,t)  =  mii(z,f)  +  mj2(z,0  (5.33) 

and  use  the  scattering  data  (5.26).  replacing  t  by  —t  in  the  second  equation  in 
(5.16)  and  adding  it  to  the  first  equation,  we  find  that 
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h{x -t)  ^  L{x ,t)  +  fh{\t-r\)L{x,r)d7  =  0  .  (5.34) 

—X 

This  result  is  knovm  as  the  Krein  integral  equation  [8], [9], [13],  and  we  have 
immediately  that  Ksi^.t)  =  L{x  ,t)  +  L(x,-t).  By  noting  that  77111(2:  .-s)  =  0  and 
that  277112(2:  .-x)  =  k{x)  yfe  find  that 

k{x)  =  —ZL{x,—x)  (5.35) 

so  that  the  local  reflection  coefficient  function,  and  therefore  the  potential,  can 
be  reconstructed  by  this  method. 

This  development  shows  that  all  the  known  solutions  of  the  inverse  scatter¬ 
ing  problem  based  on  integral  equations  can  be  related  to  the  differential 
approach  that  we  have  described  in  the  previous  sections.  The  integral  equa¬ 
tions  based  method  of  Gopinath  and  Sondhi  [5], [6]  may  be  regarded  as  using  a 
special  approach  to  the  solution  of  Krein's  equation.  This  method  is  of  impor¬ 
tance  since  the  local  impedance  is  directly  recovered  and  a  discussion  of  it  can 
be  found  in  Bruckstein  and  Kailath  [42]. 


T 
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6.  Inverse  Scattering  for  General  Media 

The  differential  inversion  methods  that  we  have  obtained  above  were  res¬ 
tricted  to  the  case  of  lossless  scattering  media.  It  is  however  possible  to  extend 
them  to  more  general  media,  where  the  wave  propagation  is  described  by  the 
two-component  system 


a 

Wjiix.t) 

dx 

Such  systems  appear  for  example  in  the  study  of  lossy  transmission  lines  and 
acoustic  media  [26], [43].  The  local  loss  function  for  this  system  is  given  by 


b(x)7iS{2r) 

Za{x) 


(6.2) 


which  shows  that  a  necessary  and  sxafficient  condition  for  losslessness  is  that,  for 
all  X. 

a(x)  =  a(x)  =  0 

l6(x)  =  ^(x)  .  (6.3) 


This  is  the  case  that  w^as  considered  in  the  previous  sections.  An  infinitesimal 
layer  of  the  scattering  medium  corresponding  to  (6.1)  is  depicted  in  Fig.  5. 

Since  the  scattering  medium  is  parametrized  by  four  different  functions 
la(  ),b  in  general  it  'will  not  be  possible  to  reconstruct  all  of  them 

from  the  pair  of  waves  and  The  reconstruction  techniques  that 

will  be  discussed  in  this  section  therefore  assume  either  that  the  parameters 
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a(2),  a{x)  and  ^(x)  can  be  expressed  as  some  function  of  5  (2),  or  that  we  have 
more  information  than  the  medium  impulse  response  at  2  =  0. 


6.1.  Reconstruction  for  a  medivun  parametrized  by  6  (2)  only 

When  the  medium  is  entirely  parametrized  in  terms  of  6(')  -  as  was  the 
case,  for  example,  for  lossless  media  -  the  layer  peeling  procedure  of  Section  3 
can  be  extended  easily  [42].  To  do  so,  note  that  when  the  medium  is  probed  by  a 
wave  with  a  leading  impulse,  the  waves  at  some  depth  2  are  of  the  form 


* 

Wj>{x  ,t)  =  (2)(5(f  — 2)  +  ^u^(2  ,b)u  (f  — 2) 

Wiix.t)  =  'Wi{x,t)u{t-x) 


(6,4) 


with 


7i?(2)  =  explfaiOd^l  .  (S.5) 

0 

In  this  case 

b(2)  =  Zy^^{x)-wi{x,x  +  )  (6.6) 

and,  since  7i?(2)  can  be  obtained  from  the  previously  reconstructed  layers, 
equation  (6.6)  may  be  used  to  compute  6(2)  for  the  next  infinitesimal  layer, 
which  in  turn  determines  0(2), a(2)  and  |3(2).  This  implies  that  (6.1)  and  (6.6) 
can  recursively  compute  the  waves  that  propagate  inside  the  medium  and  simul¬ 
taneously  recover  the  medium  parameters. 

The  requirement  that  the  medium  be  parametrized  by  &(■)  alone  might 
seem  rather  strong,  however  in  the  literature  one  often  encounters  papers  that, 
after  stating  the  problem  in  its  full  generality,  introduce  an  equivalent  assump¬ 
tion.  It  is  also  interesting  to  note  that,  in  case  the  parameters  a (2), 0(2)  and 
^{x)  depend  on  6(2)  in  a  nontrivial  way,  it  is  not  clear  how  the  integral  equa¬ 
tions  based  inversion  approaches  can  be  extended. 


w 
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6.2.  Inverse  scattering  for  a  nonsymmetric  system 

Another  example  for  which  diuerential  reconstruction  methods  can  be  dev¬ 
ised  is  when  a[x)  =  a{x)  =  0  in  (6.1).  In  this  case  the  resulting  asymmetric  two- 
component  system  is  of  the  type  considered  by  Zakharov  and  Shabat  (see.  e.g 
[16], [17]).  The  system  (6.1)  can  in  fact  always  be  reduced  to  this  particular  form 
by  performing  the  substitution 


Wi(x.t)  ■<- 

where  yniz)  is  given  by  (6.5)  and 

X 

7l{x)  =  explfaiOd^l  . 

0 

In  terms  of  these  "normalized"  variables,  (6.  i)  becomes 

7Li^) 


5 

Wjiix.t) 

dx 

b{x) 


8 

dt 


7i?(^) 

±_ 

dt 


Wp{x,t) 


(6.7) 


(6.B) 


(6.9) 


7i(«) 

which  is  now  in  the  form  of  an  asymmetric  two-component  system.  Let  us  define 
k{x)  =  and  k^{x)  =  .  (6.10) 


The  generalized  Schur  procedure  that  we  derive  next  reconstructs  the  two  func¬ 
tions  k{)  and  fc^(’)  which  are  two  independent  functions  of  the  original 
parametrization.  Thus,  unless  a[z)  =  a{x)  =  0  for  all  this  method  will  provide 
only  a  partial  reconstruction  of  the  original  medium.  Our  presentation  follows 
that  of  Yagle  and  Levy  [44]. 

In  addition  to  the  causal  pair  of  w'aves  If^(0,O  and  Wi(O.t)  that  was  used 
earlier  as  scattering  data,  it  will  be  assumed  that  we  are  also  given  a  noncavscd 


wave  pair 
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F//(0.0  =  6{t)  + 

f//(0,O  =  xu/(0.O^(-O  • 


(6.11) 


These  waves  can  be  viewed  as  obtained  by  exchanging  the  role  of  fiXO  and  Wi{-) 
and  by  reversing  time  in  a  scattering  experiment.  The  corresponding  reflection 
coeflicient  function 


R\d)  = 


is  the  (i.2)  entry  of  S  ^(o).  where 


5(«) 


Tiio) 


(6.12) 


(6.13) 


is  the  scattering  matrix  associated  with  the  medium  over  [0,“).  It  can  therefore 
be  obtained  by  probing  the  medium  from  both  ends  and  measuring  all  the 
entries  of  S  (o).  Thus,  even  though  the  knowledge  of  the  noncausal  waves 
^i?^(0.O  Is  nonphysical,  it  can  be  assumed  that  R-^{o)  or  its  anti- 
causal  inverse  Fourier  transform  R^{t)  is  obtainable.  For  the  case  of  a  lossless 
medium,  since  S{i:>)  is  unitary,  we  have 


R^{t>)  =  R{-a)  and  R^{t)  =  i?(-f)  (6.14) 


so  that  this  additional  information  is  redundant. 

The  layer-peeling  method  can  now  be  used  for  the  asymmetric  two- 
component  system  by  noting  that,  at  point  x ,  the  antic  ausal  waves  U'/'(r,f )  and 
have  the  form 


and  that 


Wl{x,t')  =  d'^x  +  t)  ■¥  •Wi{x  ,t)TJL[—X —t) 
W^{x,t)  -  Wp{x,t)u{-x-t) 


•  (6.15) 


Tnerefore,  by  using  the  system  (6.9)  ■with  the  relations  (3.3)  and  (6.16)  to  pro¬ 
pagate  both  the  caused  and  anticausal  pairs  of  waves  [Wp^x  and 

\Wi{x,t),Wp{x,t)\  simultaneously,  we  can  recover  both  A:(  )  and  k^{cdx)t)  in  a 
sequential  way. 

In  the  transform  domain  the  Riccati  equations  satisfied  by 


R{x.cS)  = 


Wi{x,u) 

Wp{x,u) 


and  R^{x,iS)  = 


Wp  (x.cj) 
fV/(x,xi) 


(6.17) 


are 


•~-/?(2,£j)  =  Zj oR{x /S)  +  k-‘^{x)R{x .cSf  ~  k{x) 


(6.18) 


~tR^{x,'S)  -=-2]  xR^  {x  ,'J)  +  k{x)R^{x  ,'Jf  -  k"{x) 
at 


These  equations  can  be  propagated  recursively  by  using  the  relations 

Urn  2j(U.ff(x,c;)  =  k{x) 


\im-Z j  uR'^{x  .cS)  =  k^{x) 


(6.19) 


(6.20) 


■which  have  the  effect  of  coupling  (6.18)  end  (6.19). 

An  integral  equations  based  solution  of  the  above  problem  can  be  found  in 


Ablowitz  and  Segur  [17]. 
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7.  Conclusionjs 

In  this  paper  we  have  obtained  differential  inversion  methods  for  identifying 
various  paremetrizations  of  lossless  and  nonlossles  scattering  media.  Crucial  in 
all  the  developments  was  the  assumption  that  the  given  scattering  data  is  noise- 
free,  therefore  the  methods  presented  are  exact  inversioTi  algorithms.  These 
methods  were  also  related  to  the  classical  approaches  of  Marchenko,  Gelfand- 
Levitan  and  Krein  which  are  based  on  solutions  of  Fredholm  integral  equations. 

Two  types  of  differential  methods  were  described.  The  layer-peeling 
method,  or  equivalently  the  Schur  recursions,  directly  exploited  the  physical 
structure  of  the  medium  to  compute  the  propagating  waves  and  to  simultane¬ 
ously  recover  its  parameters.  The  corresponding  algorithm  was  therefore  formu¬ 
lated  as  an  initial  value  problem.  A  second  set  of  differential  equations,  the 
Krein-Levinson  recursions,  were  also  derived  by  exploiting  the  structure  of  the 
Marchenko  integral  equations.  These  differential  equations  have  the  same 
dynamics  as  the  Schur  recursions,  however  they  require  certain  boundary 
values  that  have  to  be  successively  computed  by  invoking  integral  equations. 
The  algorithms  obtained  via  both  approaches  are  both  computationally  efficient 
and  numerically  stable  [40]. 

The  results  described  in  this  paper  could  be  extended  in  several  ways.  One 
of  these  would  be  their  use  for  the  propagation  of  solutions  of  certain  nonlinear 
differential  equations  by  the  inverse  scattering  transform  [l6],[36].  Also,  our 
analysis  has  been  restricted  to  physical  processes  described  by  second-order 
differential  equations.  It  would  be  interesting  to  generalize  the  differential 
approaches  discussed  in  this  paper  to  the  study  of  inverse  problems  for  more 
complex  physical  structures,  described  for  example  by  general  Hamiltonian  sys¬ 


tems. 
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FIGURE  CAPTIONS 

Fig.  1:  Discretized  wave  scattering  layer. 

Fig.  2:  Integration  path  for  the  propagation  of  the  layer-peeling  algorithm. 

Fig.  3:  Discretized  medium  associated  ■with  the  impedance  reconstruction 
procedure. 

Fig.  4:  Comparison  of  the  layer-peeling  and  layer  aggregation  methods. 

Fig.  5;  Wave  scattering  picture  for  a  general  medium. 
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